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A model inverse design problem is used to investigate the effect of flow discontinuities on 
the optimization process. The optimization involves finding the cross-sectional area distribution of 
a duct that produces velocities that closely match a targeted velocity distribution. Quasi-one- 
dimensional flow theory is used, and the target is chosen to have a shock wave in its distribution. 
The objective function which quantifies the difference between the targeted and calculated velocity 
distributions may become non-smooth due to the interaction between the shock and the discretiza- 
tion of the flow field. This paper offers two techniques to resolve the resulting problems for the 
optimization algorithms. The first, shock-fitting, involves careful integration of the objective func- 
tion through the shock wave. The second, coordinate straining with shock penalty, uses a coordi- 
nate transformation to align the calculated shock with the target and then adds a penalty propor- 
tional to the square of the distance between the shocks. The techniques are tested using several 
popular sensitivity and optimization methods, including finite-differences, and direct and adjoint 
discrete sensitivity methods Two optimization strategies, Gauss-Newton and sequential quadratic 
programming (SQP), are used to drive the objective function to a minimum. 


INTRODUCTION 

Solutions to high speed aerodynamic 
design problems generally require numerical 
solutions of the Euler or Navier Stokes equa- 
tions. These flows often contain regions of 
steep gradients such as shock waves, contact 
surfaces, and boundary and shear layers. 
Analysis of these problems are computa- 
tionally expensive. In the context of opti- 
mized design which places an even greater 
demand on computational resources, it is im- 
portant that efficient sensitivity and optimiza- 


tion algorithms are studied. In particular the 
study focuses on the effect of shock waves in 
the design optimization process. 

Recently, Frank and Shubin (Ref. 1 and 
2) studied the simple model problem of 
inviscid compressible flow through a variable 
area duct. They formulated an inverse design 
problem and investigated several design 
sensitivity techniques. Although their 
designs contained shock waves, their results 
did not indicate any adverse effects of shock 
waves on the design optimization process. 
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This appears to be due to the fact that their 
initial conditions always placed the shock at 
or very near the targeted shock position. 

The objective of the present work is to 
reexamine the model problem considered by 
Frank and Shubin and to develop efficient 
strategies for treating flows with shock 
waves. The model problem studied in Ref. 1 
and 2 offers several simplifications which are 
useful for this study. First, because of the 
quasi-one-dimensional flow approximation, 
e.g., Anderson (Ref. 3), an exact algebraic 
solution may be found, so that discretization 
errors of the numerical solution can be accu- 
rately computed. Second, accurate finite- 
volume or finite-difference solutions to the 
governing Euler equations may be found 
efficiently with shocks captured very accu- 
rately using modern computational fluid 
dynamics techniques. 

Following the inverse design problem of 
Ref. 1 and 2, we attempt to determine a 
geometry which closely approximates 
prescribed flow solutions. Inverse problems 
of this type are useful for this study, since a 
prescribed flow field can correspond to a 
known geometry, thereby giving an accurate 
measure of design errors. Specifically, the 
design variables are points defining the duct 
area distribution as a cubic spline. This 
differs slightly from Ref. 1 and 2, where 
Frank and Shubin used coefficients of B- 
sp lines. In this work B -splines are also 
implemented to compare with their results. 

Our treatment of the shock waves 
involves the use of the method of strained 


coordinates for perturbations of transonic 
flows with shock waves, introduced by 
Nixon (Ref. 4). This method has been 
applied for airfoil approximations by Stahara 
(Ref. 5) and utilized to find sensitivity 
derivatives in an airfoil optimization by Joh, 
Grossman, and Haftka (Ref. 6). 

This paper first reviews the governing 
equations for quasi-one-dimensional flow 
through a variable area duct, boundary 
conditions and the flow solutions. In the 
next sections the inverse design problem is 
formulated and the objective function is 
examined closely. The non-smooth nature of 
the objective function for flows with shock 
waves is shown to lead to convergence 
difficulties in the optimization problem. 
Methods of treating this difficulty, including 
shock-fitting and the coordinate straining are 
described. Next the optimization procedures 
used for this problem are discussed. Then 
finite-difference, direct-discrete and adjoint- 
discrete methods for computing design 
sensitivities are presented. Results and a 
discussion of the advantages and disad- 
vantages of the sensitivity and optimization 
algorithms conclude this report. 

ANALYSIS 

Governing Equations 

The governing equations for steady, 
quasi-one-dimensional flow through a duct of 
varying cross sectional area are the Euler 
equations. 


2 


puA 


'0 

(pu 2 + p)A 

► + < 

-P A x ► 

(pe 0 + p)uA 

X 

0 

Vi. 


The domain, x, varies from 0 to 1, p is the 
density, u is the velocity, A is the area, p is 
the pressure, and e Q is the total energy per 
unit mass. The equation of state for a perfect 
gas closes the system. 


and exit conditions, respectively. These 
conditions were used in Ref. 1 . 

The Euler equations are solved using 
three methods to aid in understanding the 
effects of the sharpness of the shock on the 
design process. One method is analytic and 
an exact solution is obtained. Two finite- 
volume formulations are used to get approx- 
imate solutions. The finite-volume methods 
used are Godunov and artificial viscosity. 


p = (y-l)pe, (2) 

where e is the energy per unit mass, and y is 
the ratio of specific heats assumed to be 
constant. Following the derivation in Ref. 1, 
the first and third components of (1) may be 
integrated directly, and after algebraic 
manipulation the system is reduced to a single 
ordinary differential equation in u, 

f x + g = 0, (3) 


where 


f(u) = u + 


(y-1) 2h 0 

(Y + l) u 


g(u) = 


AxY- 1 
A y + l 


(u-^H 

u 


(4) 

(5) 


Exact Solution 

The exact solution is arrived at by inte- 
grating (3) over regions where the solution is 
smooth. The result from Ref. 1 is 

Au(2h 0 - u 2 ) 1/(Y-1) = k , (6) 

where k is the constant of integration. 
Equation (6) is applied over the region from 
the inlet to the position of the shock wave, 
and again from the position of the shock 
wave to the exit. The constant, k, increases 
across the shock and is found from the 
boundary conditions at the inlet and exit for 
each application of the solution. The position 
of the shock which is the boundary of the 
right and left solution is determined by satis- 
fying the Rankine-Hugoniot relation 


and ho is the total enthalpy per unit mass. 

We specify inlet and exit velocities to get 
a unique solution to the governing differential 
equations. Velocity boundary values are 
normalized with respect to the speed of sound 
at the inlet and are 0.506 and 1 .299 for inlet 


u l u r - u * -2h °~p 

where ul and ur are the left and right values 
of the velocity at the shock. The method for 
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obtaining the exact solution is shown graphi- 
cally in figure 1 . 

Finite-Volume Solutions 

A numerical solution to (3) is obtained by 
adding an unphysical time derivative and 
marching to a steady state, 

u t + f x + g = 0. (8) 

Time integration is performed using Jameson 
four-stage Runge-Kutta. Distance along the 
duct is discretized with uj evaluated at the cell 
centers and fluxes evaluated at the cell faces. 
The spatial derivative is replaced with the 
finite-volume formulation 

u t + - fj — ~ fj - 1/2 + g j = 0. (9) 

Ax J 

Using the Godunov scheme described in Ref. 
1, the fluxes are computed according to 

fj+l, Uj, u j+1 < u* 

fj, Uj,U j+1 >U. 

f*, Uj < U* < Uj+| ’ 

max(fj,fj +1 ), Uj +1 < u* < Uj 

( 10 ) 

where fj+i = f(uj+i), etc., and * indicates 
sonic flow. This formulation yields a solu- 
tion containing a sharp shock. An alternative 
to the Godunov formulation is the artificial 
viscosity scheme used in Ref. 1 as 

f j+i/2 = j+i + fj - <*(u j+1 - Uj)] ,(11) 


f j+l/2 ~ \ 



Figure 1: Diagram of the exact solution algorithm as 
described in Ref. 1. 


where a is an artificial viscosity parameter 
which is related to the numerical dissipation. 

The exact solution resolves the shock 
precisely, while the Godunov scheme smears 
the shock over two grid points. The artificial 
viscosity parameter is a control of the amount 
of shock smearing. For example, a = 1 and 
a = 4 spread the shock over approximately 
5% and 10% of the domain respectively; a = 
12 virtually eliminates the discontinuity. The 
exact and numerical solutions for the target 
area distribution are compared in figure 2 
using a computational domain of 41 grid 
points. 



Distance Along Duct 


Figure 2: A comparison of solutions containing 

shocks. 
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DESIGN PROBLEM 

Area Description 

The design problem involves finding an 
area distribution so that the solution to (3) 
closely matches a given velocity distribution 
containing a shock. The area distribution is 
computed at discrete points in the domain 
from a set of design variables, £. Thus the 

solution to the design problem is expressed 
as a set of design variables. 

The formulation of the area distribution 
from the design variables is not unique. In 
this work, a cubic spline is fitted through the 
area at specified points along the duct as 
shown schematically in figure 3. The design 
variables to be optimized are the values of the 
area at these points. Inlet and exit areas are 
fixed at normalized values of 1.05 and 1.745. 
Clamped boundary conditions, i.e. zero slope 
at the ends are imposed to determine the 
spline uniquely. Design cases presented in 
this paper contain up to 20 design variables 
evenly distributed along the interior of the 
duct. 

An alternate method of formulating the 



Figure 3: Design variables describing the area distribu- 
tion 


area from design variables involves B- 
splines. This formulation e.g., Gerald and 
Wheatly Ref. 7, provides smooth curves and 
is implemented in this study for the purpose 
of comparing to the work of Frank and 
Shubin. 

A target velocity distribution was created 
by solving (3) using an area distribution 
described by the cubic 

A(x) = -1.39x3 + 2.085x2 + 1.05. (12) 

This area profile has the properties of A(0) = 
1.05, A(l) = 1.745, and zero slope at the 
ends. The cubic spline formulation can 
match (12) exactly and will serve as a check 
on our final design. 

A Non-smooth Objective Function 

To quantify how well a calculated veloc- 
ity distribution compares to the target, we 
define the objective function 

I(^ = ^J(u-u) 2 dx, (13) 

where u = u(x) is the target velocity distri- 
bution through the duct, and u = u(x;^) is 
the calculated velocity distribution. The 
velocities are normalized by the speed of 
sound at the inlet. In the discretization of the 
problem, the integral is approximated using 
the trapezoidal rule 
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!(£) = ( T i + r N) + 
where 


N-l 

Eii 2 I 0 4 ) 


r = (u - u)VAx , 


(15) 


and N is the number of grid points. 
Boundary conditions specify u at the inlet and 
exit to match the target exactly reducing (14) 
to 


, N-l 

I© = ^S'i 2 ' (1© 

1 i=2 

An optimum design is achieved when (16) is 
a minimum. 

Making the trapezoidal rule approxima- 
tion without regard to the position of the 
shock wave results in a non-smooth objective 
function. In figure 4, the objective function 
is drawn for a design problem involving one 
design variable and exact solutions to the 
governing equations of fluid motion. In this 
case the objective function is discontinuous. 

The jumps in the objective function result 
from the combination of the shock wave and 


the numerical evaluation of the objective 
function. The calculation of the objective 
function (16) does not involve any specific 
information about the location of the shock. 
Thus for small perturbations of the area 
distribution, provided the shock remains 
between the same grid points, the value of the 
objective function changes very little. The 
objective function is dominated by the 
differences in the target and calculated 
velocity in the segment between the shocks. 
Figure 5 is a typical plot of r = (u - u)VAx 
and demonstrates this effect. For pertur- 
bations of the area distribution which just 
moves the shock across the grid line, the 
objective function changes dramatically. On 
the other hand, a perturbation which moves 
the shock just inside the grid line will have 
almost no effect 

Without a clear picture of the nature of the 
objective function, one might try increasing 
the number of grid points to reduce round-off 
errors involved in the calculation of the 
design sensitivities. However, this will only 
create more "stairs" in the objective function 



Figure 4: Objective function plot for a one design 
variable case using an exact flow solver. 



I Distance Along Duct 

Figure 5: r and hence the objective function domi- 
nated by the region between the shock waves. 
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as now there are more cell centers for which 
the shock wave passes through. Thus 
instead of helping, the problem is worsened. 

Using a numerical solution for the flow 
solver, the resolution of the shock decreases, 
and the objective function while no longer 
discontinuous, remains non-smooth. Figure 
6 is a plot of the objective function calculated 
via the Godunov scheme, which smears the 
shock over 2 grid points. For a highly 
smeared shock, computed using the artificial 
viscosity scheme with a = 1, the objective 
function appears smooth (figure 7). The 
smooth objective function has a clear 
advantage over the discontinuous one. 



Design Variable 


Figure 6: Objective function plot for a one variable 
design problem using Godunov flow solver. 



Figure 7: Objective function plot for a one variable 
design problem using an artificial viscosity solver. 


unfortunately this comes at the expense of the 
accuracy of the flow solution. However, as 
the grid is refined, the problem will reappear. 

Shock-Fitting 

A more precise evaluation of the integral 
in (13) involves first dividing the integral at 
the location of the discontinuities and then 
applying the trapezoidal rule to each segment 
of the function. The objective function 
contains two shocks, one from the target 
distribution and one from the calculated 
distribution, thus the integral is divided into 
three segments. 



where x s and x s are the positions of the target 
and calculated shock waves respectively. 
Implied in this procedure is the knowledge of 
the precise locations of the shocks and the 
values of u and u on either side of both 
shocks. 

Applying numerical integration in this 
manner with exact solutions to (3) produces a 
well-behaved objective function (figure 8). 
(The small wiggles are due to plotting 
resolution and will be corrected in the final 
version of the paper.) 





j Design Variable j 

Figure 8: Precise numerical integration of equation 
13 yields a smooth objective function 



The distance between shocks is defined as 
Ax s = x s -x s . ( 19 ) 

The calculated velocity distribution is strained 
proportionately to the distance between 
shocks according to 


The difficulty in this method arises when 
the shock is smeared and the jump in velocity 
across the shock is not obvious. Using a 
fabricated definition for shock location may 
locate the shock in a consistent manner from 
test case to test case, but defining a right and 
left value for the velocity on either side of the 
shock is difficult. 

Coordinate Straining and Shock Penalty 

Another approach to dealing with an 
objective function which is non-smooth due 
to the presence of shock waves is the method 
of coordinate straining as developed by 
Nixon, Ref. 4. We utilize this method to 
align the target and the calculated shock 
waves, effectively ensuring the continuity of 
the objective function. 

The implementation of coordinate strain- 
ing involves defining a function, s(x), which 
equals zero at the inlet and exit, and has a 
value of 1 at the position of the target shock. 
The function is not unique, and here we 
choose from Ref. 4, 


u = u(x — sAx s ). (20) 

To apply (20) to in a discrete sense, we 
determine a grid index M such that 
X M < x i ~ sA x s < x M+1 . Then using linear 
interpolation, the strained velocity (20) 
becomes 

°i = u M + ^~' lM ( Xj-sAx,-x M ). 

( 21 ) 

Coordinate straining will transform the 
velocity distribution of figure 9 to that in 
figure 10. 



Figure 9: Typical velocity profile of target and calcu- 
lated velocity distribution. Region between shocks 
dominates evaluation of the objective function. 
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Figure 10: Strained velocity distribution. 


In the evaluation of the objective func- 
tion, using u in place of u, the dominating 
terms which exist in the region between the 
calculated and target shocks are removed. In 
the process of removing the staircase, the 
objective function becomes very flat near the 
minimum, thus slowing convergence (figure 
11). Using this technique, we rely complete- 
ly on the small differences in velocities 
outside the shock region to drive the area to 
the target. Often this is enough to improve 
the design, but not enough to achieve the best 
possible one. 

To shape the objective function to capture 



Figure 11: Strained objective function. 



Design Variable 

Figure 12: Modified objective function, o= 5. 

the valley of figure (8), a shock penalty 
proportional to the square of the difference 
between the calculated and target shock wave 
is added to the strained objective function, 
yielding, 

, N-l 

I = -[X f i Z + a ^s“ X s) 2 ]> ( 22 ) 

L i=2 

where q = (u,- Uj)VAx, and a is a positive 
constant. Values of a can be chosen so that 
in the first design iteration (22) equals (13). 
Figure 12 shows the objective function 
modified by coordinate straining and shock 
penalty. 

OPTIMIZATION ALGORITHMS 

Two optimization strategies were used to 
solve the design problem. One is the Gauss- 
Newton method which requires some deriva- 
tion for application to the objective function 
defined by the coordinate straining and shock 
penalty technique. The second, an SQP 
algorithm, can be applied without any special 
treatment to both the objective functions. 
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Gauss-Newton Method 

In this section the Gauss-Newton opti- 
mization algorithm is applied to the function 
defined by (22). A necessary condition for a 
minimum is 


ai 


N-l 


= S ? i^- a ^ s - x s)lr =0 ’ (23) 

i=2 


for j = 1,.., n where n is the number of 
design variables. The solution to the opti- 
mization problem requires finding the root 
that satisfies VI(^) = 0. Given the solution at 
the £th iteration, I; , we can find 
= + by expanding VI(^ +1 ) in a 

Taylor series and retaining only the first term. 
Thus using (23) we proceed via a Newton 
method to get 


comprises the Gauss-Newton method. The 
system is then 


SA5 t [xJHib <#!?]= 

k=l i=2^k 


5fj dx s 
i=2 


(25) 


To implement (25) we must accurately 
compute drj/d^j. In the following sections 
we investigate several methods to compute 
these derivatives. 


SQP Method 

The sequential quadratic program method is 
described in detail in many references, e.g. 
Haftka and Giirdal, Ref. 8, and no details are 
provided here. The algorithm used in this 
work is developed by Schittkowski, Ref 9. 


N-l 


^ ?i §' _CT(Xs_Xs) ^F i + 

i=2 


n - N ~ L a?i dx- 


+ r ; - 


a 2 h 


•i 


dx. dx. . d 2 x< 

= 0 . 


■te 

(24) 


SENSITIVITY METHODS 

Finite-Differences 

The optimization routines require either 
d rj/d^j and dxjd^j or dl/d^j which can be 
computed using finite-differences. In this 
work, first order forward and backward 
approximations. 


Applying (24) to each of n design variables 
results in a linear system of n equations 
which drives the design variables to their 
optimal value. Near the minimum, q and 
(x s — x s ) are small, and the second derivative 
terms may be neglected. This avoids the 
computation of the second derivatives and 






(26) 


dsj 

~ r i(^l*^2»”‘»^j — A^j,*",^ n ), (27) 
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and second order central difference approxi- 
mations, 

^ = ^[r i < 5 1 .5 2 .-4j+A5 j ,".,5 n ) 

-nfti.fe-.4j-A5j.---W. (28) 

were investigated. Numerical results did not 
indicate a clear advantage to any of the three. 
The central difference approximation requires 
2N additional forward solutions, whereas the 
one sided differences require only N addi- 
tional forward solutions. Results presented 
in this paper use the forward difference 
approximation. 

We concluded from a study that a compu- 
tational domain containing at least 41 grid 
points and Ai; = 1CH are sufficient for 
computing the design sensitivities accurately. 

Direct Discrete 

This method computes the design sensi- 
tivities by applying the chain rule of differen- 
tiation to the discrete governing equations. 
While this method is cheaper than the finite- 
difference approach, it is more involved to 
implement. Further, the calculation of the 
sensitivities requires knowledge of the flow 
solving algorithm. 

In this formulation, we distinguish 
between the flow variables, u, and the design 
variables The flow variables are the values 
of the velocity at the grid points in the 
domain, and are themselves functions of the 
design variables. In addition to u, we 


include the shock position, x s , which is also 
dependent on Considering the strained 
and shock penalized objective function we 
have, 

I = I[u($),x s (5)]. (29) 


Applying the chain rule to (29) to compute 
the sensitivity, we have 


ai ai du ai dx s 
34, au^ j + 3x,3^ J ' 


(30) 


where 


ai_ 

du 


_ai_ _ai_ 

dtij du 2 


91 

du 


N. 


(31) 


and 


du 


duj du 2 


du 


-iT 


N 


j 


(32) 


The straining function is defined in (21) 
and in general is a function of u and x s , 

u = u[u(£),x,(€)]. (33) 


Differentiation with respect to the jth design 
variable yields 


du _ du du du dx s 

3|; ‘3^ + 3^’ 


( 34 ) 


where 


11 



duj 

du t 

du! 

du 2 

duj 

du 2 

du t 

du 2 

du* 

duN 

duj 

du 2 


du _ du t du 2 du N 

*j 1*I Wi aI7 ’ 


x s [u(S)] = X js + (u. - u is ) . 

U js + 1-Ujs J 


The derivative of the shock position with 
respect to the jth design variable is thus 


dx s _ dx s du 
du d^’ 


where 


afl Jaa, au 2 ao N 


dx s dx s 


The shock position is defined where the 
Rankine-Hugoniot relation (7) is satisfied. In 
a numerical solution where ul and ur are not 
clearly defined, we take the position of the 
shock to be the point in a steep compression 
where the Mach number is one, which may 
be written as 


u(x s ;£) = u* = 2h 0 ^— \ 

Y Y + l 


Discretized, we locate the shock by linear 
interpolation 


u s = ujs + - 


u js+l~Ujs 


(Xg — Xj s ) = u* , (39) 


where the shock lies between the js and js+1 
grid points. Solving for the shock position 
yields 


ax i _fax s dx s 


du a Ul du. 


au N J' 


Substituting (34) and (41) into (30) gives us 
an expression for the design sensitivities 


ai _ V T du 
^i = 


where 


yT^f— +— ^]+— ^M44) 

du^du dx s du J dx s du 


From equations (21), (22), and (40) we 
can obtain analytic expression to evaluate all 
the derivatives in (43) with the exception of 
du/d^j. The direct discrete method applies 
the chain rule to the governing equations 
discretized by some numerical scheme. We 
have N discretized equations that are, in 
general, functions of the flow variables and 
the design variables. These equations are 
given the symbol w. 
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w i = w 1 (u 1 ,.. M u N ^ 1 ,...4 n ) = 0 
w 2 = w 2 (u 1 ,.. M u N ^ 1 ,..4 n ) = 0 

w N=w N (u 1 ,.. M u N ,^ 1 ,...^ n ) = 0. (45) 

Differentiating (45) with respect to the jth 
design variable, we find 


Equation (50) must be solved n times to 
find the sensitivities with respect to all the 
design variables. For other problems where 
N is extremely large, it may not be possible 
to store the L-U decomposition of J. In such 
a case solving (50) n times would be 
expensive. 


dw | dw ] , 0u 

*j"UJL *r 0 ' 


m 


where 


w =K w 2 ... w n j t . 


(47) 


dw 


dw : dw 2 dw N 

% Hi % 


iT 


, (48) 


and 


3wj 

dw ( 

9w, 

duj 

3u 2 

du N 

dw 2 

dw 2 


du l 

du 2 

... 

9u N 

dw N 

3 wn 


dui 

du 2 

... 

du N 


j-i^i = 

ou 


(49) 

We can find 3u/3^j by solving the linear 
system 



(50) 


Adjoint Discrete 

A third method to compute design sensi- 
tivities from the discrete equations is the 
adjoint method. We define an augmented 
objective function 


I -I + A, t w , ( 51 ) 

where I is defined in (22), AT is a row vector 
of Lagrange multipliers and w is the column 
vector of discretized governing equations. 
The sensitivity with respect to the jth design 
variable is 

dl _ dl T dw 

wrwi +x w; (52) 


The first term on the right hand side has 
been expanded in the derivation of the direct 
discrete method (equation (43)). The deriva- 
tive of the governing equations has also been 
expanded in equation (46). Substituting (43) 
and (46) into (52) yields 


ai , _ 

d $i 


-ti 


r dw^ 




• (53) 
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We rearrange (53) to collect terms multiply- 
ing du/^j to get 


af / 



aT • (54) 


The vector of Lagrange multipliers is arbi- 
trary. If A. is chosen such that 


J T X = -v, 

then the sensitivity is 
dl _ 91 _ . t J 9w | 

I* f: 


(55) 


(56) 


The advantage this method has over the 
direct discrete method is evident with very 
large systems. Equation (55) has to be 
solved only once to calculate all the Lagrange 
multipliers. The sensitivities maybe then 
computed inexpensively using (56). 


Comparison 

For a design case involving three design 
variables and the artificial viscosity flow 
solver, sensitivities have been recorded in 
tables 1, 2, and 3. Each table contains sensi- 
tivities calculated by each method. Table 1 
corresponds to the discontinuous objective 
function, table 2 uses coordinate straining, 
and table 3 uses coordinate straining with a 
shock penalty. We see that all three methods 
produce similar results for the design 
sensitivities. 


Finite- 
Difference 


Direct 

Discrete 


Adjoint 

Discrete 


ai/dS, 

ai/^2 

ai/a^ 3 


-4.914e-2 

6.200e-2 

-2.995e-l 


-4.921e-2 
6.23 le-2 
-3.003e-l 


-4.92 le-2 
6.23 le-2 
-3.003e-l 


Table 1: Sensitivity comparison using the discontin- 
uous objective function. 


Finite- Direct 

— Di fference Discrete 

dl/a^, 2.312e-3 2.305e-3 

ai /^2 4.012e-4 3.996e-4 

-2.514e-3 -2.534e-3 


Adjoint 

Discrete 

2.305e-3 

3.996e-4 

-2.534e-3 


TabIe 2: Sensitivity comparison using the coordinate 
strained objective function. 


Finite- 

Difference 


Direct 

Discrete 


Adjoint 

Discrete 


ai/5^1 -1.355e-2 -1.369e-2 -1.369e-2 

dl/d% 2 -5.326e-2 -5.260e-2 -5.260e-2 

__ _jJ41e-l -8.741e-l 

Table 3: Sensitivity comparison using the coordinate 
stramed and shock penalty objective function. 


RESULTS 

Design cases were run varying the 
number of design variables, the flow solution 
algorithm, the method of calculating the 
sensitivities, and the optimization routine. 
The initial values of the design variables 
describe the area distribution shown in figure 
13. The initial distribution places the shock 
wave significantly far away from the target 
shock so that the interaction of the flow 
discontinuity and the discretization of the 
flow field on the optimization process 
becomes important. Figure 14, 15, and 16 
show the initial velocity solution in compari- 
son to the target for each of the flow solving 
algorithms. Average CPU times on a SGI 
Iris 340 VGX to compute the velocity field 
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Figure 13: Initial area Distribution 



Figure 14: Initial velocity profile using exact flow 
solver. 



Figure 15: Initial velocity profile using Godunov 
flow solver. 


i.e., solution to (3) are listed in table 4. The 
exact solution was computed in a very small 
amount of CPU time, less then 0.5 sec. 



Figure 16: Initial velocity profile using artificial 
viscosity flow solver. 



Ave. CPU 


(sec.) per sol'n 

Exact Solution 

0 

Godunov Solution 

8 

Artificial Viscosity Solution 

3 


Table 4: CPU time for flow field solutions. 


Shock-Fitting 

The shock-fitting objective function was 
tested using the finite-difference method with 
the SQP optimization routine. As mentioned 
earlier, the objective function is tailored for 
the case where the jump in velocity across the 
discontinuity is clearly defined. For this 
reason only exact flow solution results are 
presented. In comparison, the design process 
is attempted using the discontinuous objective 
function. Tables 5, and 6 list the initial, 
final, and targeted values of the design 
variables for optimization of the discon- 
tinuous and shock-fitted objective functions 
using 3 and 7 design variables. We see that 
the shock-fitted objective function will 
converge to the target area distribution 
whereas the discontinuous objective function 
converges to the incorrect solution. The SQP 
optimization algorithm allows the user to 
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specify maximum and minimum limits on the 
design varaibles. In the case of the 
discontinuous objective funciton, many 
design variables were driven to these limits. 


^initial 

Discontin 

-uous 

^final 

Shock- 

fitted 

^final 

^target 

i.i 

1.1110 

1.1822 

1.1586 

1.2 

1.1824 

1.4010 

1.3975 

1.4 

1.3830 

1.7662 

1.6364 


Table 5: Shock-fitted results for 3 design variables 


^initial 

Discont. 

£final 

Shock- 
fitted 
£ final 

^tareet 

1.0643 

1.0879 

1.0718 

1.0799 

1.1000 

1.1404 

1.1644 

1.1586 

1.1471 

1.3000 

1.2625 

1.2700 

1.2000 

1.2000 

1.3988 

1.3975 

1.2658 

1.2500 

1.5164 

1.5251 

1.4000 

1.4000 

1.5958 

1.6364 

1.6202 

1.8000 

1.7139 

1.7151 

Table 6: Shock-fitted results for 7 design variables 


Comparison of Sensitivity Algorithms 

In the following paragraphs the design 
processes using the finite-difference, direct 
discrete, and adjoint discrete methods for 
computing design sensitivities are compared. 
Specifically, the cost measured in CPU time 
for computing the sensitivities and time to 
compute the design are used as a basis for 
comparison. In this comparison, the artificial 


viscosity algorithm is used for the computa- 
tion of the flow solution and the Jacobian, J. 

Table 7 lists the pertinent data for this 
comparison. In table 7, F-eval refers to the 
number of flow field solutions needed to 
converge the design. This number is equal to 
one plus the number of forward evaluations 
needed per iteration since an extra solution is 
required to check for convergence. 
Sensitivity CPU refers to the average time to 
compute one design sensitivity dl/d^. CPU 
time indicates the total time required to 
complete the optimization. 

As expected the direct discrete method is 
much more efficient than the finite difference 
method. 

CONCLUSIONS 

In this paper we have developed several 
algorithms to calculate design sensitivities. 
In the application of these methods, compli- 
cations involving the interaction between 
flow discontinuities and the discretization of 
the flow field cause the design process to fail. 
We have shown that both shock fitting and 
the coordinate straining and shock penalty 
technique can overcome these difficulties. 



No, of Design 

No. of 

■n 

Sensitivity 

CPU 


Variables 

Iterations 


CPU (sec) 

Time 

Finite-Difference 

3 

13 

53 

2.79 

179 

Finite-Difference 

7 


241 

3.54 

937 


19 

39 

1561 

3.28 

2323 

Direct Discrete 

3 

14 

15 

.010 

97 

Direct Discrete 

7 

32 

33 

.0057 

229 


19 

39 

40 

.0042 

243 

Adjoint Discrete 

3 

14 

15 

.010 

98 

Adjoint Discrete 

7 

32 

33 

.0043 

228 

Adjoint Discrete 

19 

39 

40 

.0021 

241 


Table 7: Comparison of Sensitivity Algorithms using the artificial viscosity solver 
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Sensitivity methods are easily adaptable to 
coordinate straining with a shock penalty 
which makes the optimization algorithm 
efficient and robust for cases involving flow 
discontinuities. The finite-difference algo- 
rithm, direct discrete and adjoint discrete 
method provide similar design results. The 
finite-difference method solves the forward 
problem N+l times per iteration and is 
impractical where the flow calculation is 
expensive. The direct methods, though more 
involved to implement, have significant time 
savings. 
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